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1.  INTRODUCTION 

1.1  RAY  TRACING  & MINIMUM  GROUP  PATHS: 

One  of  the  standard  techniques  for  the  study  of  electro- 
magnetic wave  propagation  in  inhomogeneous  media  is  ray 
tracing,  which  is  based  upon  the  geometrical  optics  formula- 
tion. (For  reviews  of  ray  theory  especially  applied  to  the 
ionosphere,  see  Kelso,  1964;  Yeh  and  Liu,  1972).  This  tech- 
nique allows  us  to  treat  the  wave  propagation  problem  as 
a problem  in  mechanics  in  which  our  objective  is,  essentially, 
to  calculate  the  trajectory  of  a fictitious  point  (or  particle) 
that  travels  in  space  with  a velocity  equal  the  local  group 
velocity  of  the  electromagnetic  wave.  The  knowledge  of  the 
ray  trajectory  also  enables  us  to  calculate  other  parameters 
(such  as  group  path,  phase  path,  and  qround  range)  which 
may  or  may  not  be  relevant  to  a particular  problem  in 
electromagnetic  wave  propacjation . 

While  the  form  of  the  equations  which  govern  the  ray 
trajectory  (i.e.  the  ray  equations)  are  fairly  general 
most  of  the  ionospheric  applications  of  the  ray  equations 
are  based  upon  a spherical  co-ordinate  system  in  which  case 
the  ray  equations  are  usually  written  _n  a form  referred  to 
as  Haselegrove ' s equations  (Haselegrove , 1957;  Yeh  and  Liu, 
1972).  In  the  present  work,  we  have  restricted  our  attention 
to  a two  dimensional  case  in  which  the  rays  are  assumed  to 
propagate  in  one  member  of  a family  of  vertical  planes  con- 
taining the  straight  azimuthal  lines  directed  outward  from 
the  source.  Furthermore  it  is  assumed  that  the  medium  in 
which  the  rays  propagate  (that  is,  the  ionosphere)  has 
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properties  which  vary  in  only  the  radial  and  angular  co- 
ordinates (see  Figure  1).  In  other  words,  the  ray  propaga- 
tion equations  become  two  dimensional.  It  is  furthermore 
assumed  that  the  geomagnetic  field  may  be  neglected  in 
which  case  the  dispersion  surface  of  the  rays  becomes 
isotropic  (or,  equivalently  the  wave-normal  and  group 
velocity  vectors  are  parallel) . With  these  conditions  in  mind, 
the  ray  equations  may  be  written  in  the  form  of  a set  of 
coupled  ordinary  first  order  differential  equations  (i.e. 

Yeh  and  Liu,  1972) ? 


(1) 


de  „ ^e_ 

dT  rn2 

dqr  _ 1 9n  + o d0  _ I9n  + °6 
dx  n 9r  0 dx  n9r  rn2 


(2) 

(3) 


dq9  _ 1_  9n  _ ^ dr  _ i 9n  _ q0qr 
dx  rn  90  r dx  rn  96  rn2 


with 


n = phase  refractive  index 


o = nk  • r 
r 

A A 

a q = nk*0 


and  with  dx  representing  the  differential  element  of 
phase  path.  Two  additional  equations  have  also  been 
used  in  order  to  allow  us  to  keep  track  of  the  cumulative 
group  and  phase  paths  as  the  ray  progresses.  These  additional 


w 


equations  are: 


" = 1 
dx 


and 


(5) 


dP ' 
dx 


(6) 


with  the  first  equation  representing  the  phase  path  and 
the  second  equation  representing  the  group  path.  The  set 
of  initial  conditions  associated  with  equations  ( i through 
(6)  is  based  upon  the  assumption  that  there  is  a free  space 
region  (n=l)  between  the  surface  of  the  Earth  (source  loca- 
tion) and  the  base  of  the  ionosphere.  Thus,  the  initial 
conditions  are: 


r = 6370  km  = radius  of  Earth 

0 =0°  (selection  of  source  location  as 

angular  reference) 

o = sin^  0 ( 4>  0 = initial  ray  elevation  angle) 

°0  = cos 
and,  of  course 

P = P'  = 0 


(7) 

(8) 

(9) 

(10) 

(11) 


The  net  result  of  using  these  initial  conditions  in 
connection  with  a given  model  for  the  phase  refractive 
index,  n,  is  that  the  only  remaining  degree  of  freedom  for 
a ray  of  some  specified  frequency  is  provided  by  the  initial 
elevation  anqle  of  the  ray  (<{>„)  and  consequently  the  value 
of  group  path,  P',  which  is  computed  when  the  ray  has  left 
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the  source,  reflected  from  the  ionosphere,  and  returned 
back  to  the  surface  of  the  Earth  is  determined  by  <f0*  Of 
course,  certain  choices  of  tp0  may  result  in  the  ray  pene- 
trating the  ionosphere.  Nevertheless,  for  those  values 
of  <j>0  which  result  in  rays  returning  to  the  surface  of  the 
Earth  there  will  be  some  value  of  <)>„  which  results  in  a 
minimum  group  path.  The  physical  significance  of  minimum 
group  path  rays  is  that  they  constitute  the  leading  edge 
of  backscatter  ionograms  (as  discussed  by  Rao  et  al . , 1976) . 
This  particular  curve,  in  turn,  forms  the  basis  for  the 
selection  of  the  ionosphere  model  parameters. 


1.2  LQP  MODEL  & IONOSTATES 


In  earlier  work  (Rao  el;  al . , 1976)  we  used  the  quasi- 
parabolic model  to  represent  the  ionospheric  electron 
density  profile  (Croft  and  Hoogasian,  1968)  in  accordance 
with : 


NeU)  - »B[1 


2 r 2 

(1_!E)  (_b)  j r<r<r(-r-b— ) 

t,  r b m 


r m 


rb-y: 


m 


where 


(12) 


N 

m 


rb 

Ne(r) 


maximum  electron  density 
radius  to  peak  of  ionosphere 
layer  semi-thickness  = rm-rb 
base  radius 

electron  density  at  a geocentric  distance  of  r 
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Equation  (12)  may  also,  however,  be  used  in  order  to  obtain 
an  expression  for  the  phase  refractive  index,  n,  since 


n 


2 _ 


= 1 - 


80. 6N, 


(13) 


In  this  case,  it  becomes  more  convenient  to  use  the  critical 
frequency  (fc)  parameter  instead  of  the  maximum  electron 
density  parameter  (Nm) . These  parameters  are  related,  in 
the  international  system  of  units,  by 


fc2  = 80.6Nm 


(14) 


Thus,  the  refractive  index  at  some  given  geocentric  dis- 
tance, r,  may  be  computed  by  specifying  the  three  parameters 

rb ' rm  and  fc- 

ln  the  present  work,  however,  we  have  attempted  to 
extend  the  quasi-parabolic  model  to  allow  the  values  of 
rb'  rm'  and  fc  to  var¥  with  the  angular  co-ordinate,  0. 
Formally,  the  electron  density  may  then  be  written  as 

r-r_(0>  2 r.(9)  2 

Ne  ( r , 0 ) = Nm  ( 0 ) { 1 - [ — — ] [_ ] } 

Ym ( 9 ) r 

(15) 


which  introduces  a corresponding  angular  dependence  in  the 
phase  refractive  index.  Such  an  angular  dependence  is 
desirable  in  situations,  for  example,  of  simulating  ion- 
ospheres in  the  trough  region.  Nevertheless,  for  a fixed 
value  of  0,  equation  (15)  has  the  form  of  an  ordinary  quasi- 
parabolic model  of  electron  density  and  we  have  accordingly 
decided  to  designate  this  model  as  the  LQP  (or  locally 
quasi-parabolic)  model. 
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In  most  of  the  cases  we  will  present  it  will  be 

assumed  that  the  angular  dependence  of  the  model  parameters 

consists  of  a nominal  value  which  is  given  by  the  value  of 

the  parameter  directly  above  the  source  (9=0)  plus  a constant 

gradient  in  the  0 direction  multiplied  by  the  angular  dis- 

drh 

placement  away  from  the  source  (e.g.  rb  (e)  = + dT  ' 9 = 

rbo  + GRB • 0 ) . Since  there  are  three  parameters,  rbo,  rmo, 
and  fco»  as  well  as  three  gradients  (one  in  each  parameter) , 
there  is  a total  of  six  parameters  needed  to  describe  the 
LQP  model.  These  six  parameters  may  be  considered  to  define 
a basis  for  a six  dimensional  space  in  which  a simultaneous 
set  of  the  six  parameters  may  be  thouqht  of  as  constituting 
a vector  (the  ionostate)  which  we  may  define  as: 

t = [RB,  RM,  FC,  GRB,  GRM,  GFC ] (16) 

with 


RB 

— 

rbo 

RM 

= 

rmo 

FC 

= 

fco 

GRB 

= 

drb/d0 

GRM 

= 

drm/d0 

GFC 

= 

dfc/d0 

1.3  METHOD  OF  SOLUTION 

Having  defined  the  ionostate  above  (equation  16)  the 
problem,  stated  in  its  simplest  terms,  becomes  now  to  deter- 
mine the  six  components  of  the  ionostate  in  some  optimum 
sense.  The  criterion  which  we  have  selected  is  that  the 
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backscatter  leading  edge  which  would  be  observed  if  the 
ionosphere  were  LQP  should  match,  as  closely  as  possible, 
the  backscatter  leading  edge  as  experimentally  measured 
by  the  ionospheric  sounder  operating  at  the  same  set  of 
frequencies  and  observing  the  actual  ionosphere.  More 
precisely,  if  we  denote  the  set  of  sounding  frequencies 
by  f^  (for  i=l  to  N)  and  denote  the  measured  backscatter 
leading  edge  by  P'  (f^) , while  denoting  the  LQP  model  back- 
scatter  leading  edge  as  P'(£,f^),  the  problem  may  be  stated 
as : 

Find  £ such  that 

N 

e2  (t)  = ^ l [P ' (f^-P'  (f  ,f±)  ]2  (17) 

i=l 

is  a minimum.  A discussion  of  the  mechanics  of  computing 
the  appropriate  ionostate,  f,  will  be  deferred  until 
section  3. 
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2 . RESULTS 

2.1  TWO  VARIABLE  CASE 

While  the  ionostate  as  previously  defined  has  six 
components,  our  present  computer  program  has  the  capability 
of  distinguishing  between  "fixed"  and  "variable"  ionostate 
components.  A fixed  ionostate  component  implies  that  the 
reduction  in  the  mean  squared  error  (e2)  as  defined  by 
equation  (17)  may  not  be  accomplished  through  any  change 
in  this  particular  component.  Instead,  the  program  may 
only  reduce  the  mean  squared  error  by  adjusting  the  vari- 
able ionostate  components.  In  this  section  we  will  be 
considering  two  simultaneous  "variable"  components. 

Furthermore,  in  order  to  provide  a reference  by  which 
to  determine  the  efficacy  of  the  program  it  is  necessary 
to  simulate  the  measurement  of  the  backscatter  leading 
edge  in  an  ionosphere  of  known  refractive  index.  For 
this  purpose,  we  assumed  that  the  ionosphere  was  locally 
quasi-parabolic  and  was  completely  specified  by  the  iono- 
state,  £„  where 

f0  = [RB;  RM;  FC,  GRB ; GRM;  GFC]  = 

[6570  km;  6720  km;  5 MHz;  0 km/rad;  0 km/rad;  2MHz/rad.] 

(18) 

-V 

Based  upon  this  we  then  proceeded  to  calculate  the 

minimum  group  paths,  P'0  (f ^)  , corresponding  to  a backscatter 

ionogram  at  discrete  sounding  frequencies.  We  then  "esti- 

mated"  an  initial  ionostate  E,  ^ which  was  found  to  result  in 

-*■ 

minimum  group  paths,  P*  ( f i ) . The  estimation  of  £ , however. 
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is  not  always  a simple  matter  and  the  reader  may  wish  to 
refer  to  section  4 for  a more  thorough  discussion.  Using 
the  mean  squared  error  e2 , defined  in  equation  (17)  as  a 
criterion  of  optimality,  we  attempted  to  reconstruct  the 
original  ionostate.  £„• 

As  an  example,  the  cases  illustrated  in  Figure  2 result 
from  specifying  the  RM,  FC,  GRB,  and  GRM  components  of  Sj 
as  "fixed",  thereby  reducing  the  problem  to  a two  dimensional 
minimization.  In  all  these  cases  the  fixed  components 
were  set  to  their  correct  values  (i.e.  the  values  of  the 
corresponding  components  of  !0).  Using  three  different 
sets  of  starting  values  (marked  by  x's  in  Figure  2)  for 
the  remaining  components  (RB  and  GFC) , the  program  was 
allowed  to  continue  in  an  iterative  manner  until  the  mean 
squared  error  became  less  than  1 km2 . The  underlined 
numbers  in  this  figure  represent  the  mean  squared  error, 
and  the  dotted  curve  corresponds  roughly  to  the  contour 
curve  along  which  the  mean  squared  error  is  approximately 
7.  A somewhat  more  complete  set  of  contour  curves  is  shown 
in  Figure  3.  As  can  be  seen  in  comparing  Figures  2 and  3, 
the  trajectory  of  the  point  (RB,  GFC)  proceeds  primarily 
along  the  direction  of  a narrow  valley  in  the  mean  squared 
error  surface. 

Figure  4 shows  a similar  convergence  when  RM  & GFC  are 
used  as  variable  ionostate  components. 

2.2  SIMULTANEOUS  DETERMINATION  OF  3 GRADIENTS 

In  this  section,  we  assume  that  the  overhead  electron 
density  profile  is  known.  Therefore,  the  values  of  RB,  RM, 
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and  FC  may  be  determined  by  finding  the  quasi-parabolic 
model  which  best  fits  the  overhead  profile.  The  inversion/ 
synthesis  program  may  then  be  used  to  determine  the  gradient 
quantities  in  the  LQP  model  (i.e.  GRB , GRM,  and  GFC) . In 
order  to  simulate  this  use  of  the  program,  we  started  out 
by  assuming  the  ionosphere  to  be  described  by  an  LQP  model 
with  the  following  ionostate: 

RB=6570  km;  RM=6720  km;  FC=5  MHz;  GRB=25  km/rad;  GRM=20  km/rad; 
and  GFC=2  MHz/rad. 

Schematically,  this  particular  model  would  look  somewhat 
like  Figure  5.  Thus,  the  gradient  in  RB  caused  the  base 
of  the  ionosphere  to  recede  from  the  surface  of  the  Earth 
as  distance  from  the  transmitter  increased.  Similarly, 
the  gradient  in  RJ1  caused  the  peak  height  to  recede  from 
the  surface  of  the  Earth  (although  at  a rate  somewhat  less 
than  that  involving  RB) . Thus,  the  ionosphere  became 
narrower  with  distance  away  from  the  transmitter  site.  The 
gradient  in  critical  frequency  (GFC) , however,  tended  to 
make  the  electron  density  greater  with  increasing  distance 
from  the  transmitter.  As  before,  the  first  step  was  the 
generation  of  a backscatter  leading  edge  in  an  ionosphere 
with  an  ionostate  of  {6570,  6720,  5,  25,  20,  2}.  The  re- 
sults were  computed  for  sounding  frequencies  of  6.5,  7.0, 
and  7.5  MHz,  and  are  summarized  below: 


FREQUENCY  (MHz) 

MINIMUM  GROUP  PATH  (KM) 

6.5 

1081.2 

7.0 

1172.3 

7.5 

1265.4 

INCREASING 

ELECTRON 

DENSITY 


TRANSMITTER 


EARTH  SURFACE 


CENTER  OF  EARTH 


Figure  5.  Combined  effects  of  three  simultaneous  gradients 
upon  ionosphere  model. 
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The  inversion/synthesis  program  was  then  run  in  the 
inversion  mode  for  two  separate  cases.  In  the  first  case 
an  initial  ionostate  of  {6570,  6720,  5,  15,  10,  1.9}  was 
selected  with  the  last  three  components  specified  as  variable 
ionostate  components.  The  mean  squared  error  resulting 
from  this  particular  ionostate  was  found  to  be  rather  small 
(3  km2)  and  thus  the  ray  tracing  was  performed  with  a fairly 
small  step  size  (.5  km).  The  following  figure  (Figure  6) 
illustrates  the  convergence  of  the  ionostate  in  terms  of 
the  mean  squared  error.  The  top  figure  (6a)  may  be  thought 
of  as  the  projection  of  the  three  dimensional  ionostate 
trajectory  upon  the  GRB-GFC  plane,  while  the  lower  figure 
(6b)  may  be  thouqht  of  as  the  projection  of  the  same  tra- 
jectory upon  the  GRM-GFC  plane.  The  small  underlined 
numbers  beside  the  trajectory  represent  the  mean  squared 
error  which,  in  this  case,  was  reduced  from  3 km2  to  2 km2 
after  the  first  iteration  and  from  2 km2  to  1 km2  after 
the  second  iteration  (all  mean  squared  errors  have  been 
rounded  to  one  significant  digit).  Based  upon  past  experience, 
it  would  seem  that  the  mean  squared  error  may  be  reduced 
beyond  1 km2 , but  this  would  require  a significantly  re- 
duced step  size  in  the  ray  tracing  subroutine. 

A somewhat  more  convincing  demonstration  of  the 
simultaneous  determination  of  the  three  gradients  is 
provided  by  a second  case.  In  this  instance,  we  selected 
an  initial  ionostate  of  {6570,  6720,  5,  40,  30,  1.8},  and 
found  an  initial  mean  squared  error  of  43  km2.  As  a result 
of  the  relatively  large  mean  squared  error  (when  compared 
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to  the  previous  case)  we  increased  the  step  size  from  .5  km 
in  the  previous  case  to  2 km.  The  ionostate  trajectory 
for  this  case  is  displayed  in  Figure  7 using  the  same  con- 
vention as  were  used  in  the  previous  case. 

In  conclusion,  the  study  of  the  simultaneous  inversion 
of  three  gradients  has  shown  the  applicability  of  the 
present  program.  In  general,  however,  these  inversion 
runs  have  been  rather  expensive  in  terms  of  computer  time. 
Efficient  use  of  the  program  for  this  type  of  problem 
would  then  seem  to  require  judgment  regarding  the  step 
size  in  the  ray  tracing  program. 

2.3  THE  EFFECT  OF  UNDERLYING  IONIZATION 

Several  computer  simulations  were  run  in  order  to 
demonstrate  the  influence  of  underlying  ionization  upon 
both  the  generation  of  the  backscatter  ionogram  and  the 
profile  resulting  from  the  ionogram  inversion.  These  simu- 
lations required  a few  minor  modifications  in  the  ray 
tracing  subroutine  which  appears  in  the  inversion/synthesis 
program.  The  net  effect  of  these  modifications  was  to 
introduce  a quasi-parabolic  electron  density  model  (with 
a base  height  of  130  km,. and  a semithickness  of  100  km)  which 
would  combine  additively  with  the  main  LQP  model  (6570,  6720, 
5,  0,  0,  2).  As  an  example  of  the  overall  electron  density 
model  which  results  when  the  critical  frequency  of  the 
underlying  ionization  is  set  at  .75  MHz,  the  following 
graph  (Figure  8)  may  be  considered.  This  graph,  of  course, 
corresponds  to  the  profile  which  would  exist  above  the 
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35  GRM  (km/rod) 
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Figure  7.  Three  component  (GRB  ,GRM,  S.GFC)  convergence  of 
MSE  with  a large  initial  error. 


o 

o 

O 

o 

o 

O 

N- 

ID 

IT) 

CO 

ID 

<D 

(W>1)  30IWlSia  0IH1N30039 

Fiqure  8.  Combined  electron  density  model  (with  underlying  layer). 
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sounder.  This  distinction  is  necessary  since  the  main 
LQP  model  contains  a gradient  in  critical  frequency,  while 
the  underlying  layer  is,  in  essence,  a spherically  symmetric 
electron  density  model. 

In  order  to  determine  the  effects  of  the  underlying 
layer  upon  the  present  problem,  we  compared  the  results 
of  simulated  ionogram  generation  (synthesis)  in  both  the 
presence  and  absence  of  the  underlying  layer.  The  minimum 
group  paths  as  measured  at  6.5,  7.0,  and  7.5  MHz,  in  the 
absence  of  the  underlying  layer,  provided  a suitable 
reference,  and,  in  fact,  the  values  of  minimum  group  path 
as  measured  at  these  three  frequencies  were  1076.3  km, 

1165.4  km  and  1256.3  km,  respectively.  Using  the  previously 
described  model  for  the  underlying  ionization,  various 
values  of  the  critical  frequency  of  this  layer  were  used 
in  order  to  simulate  various  degrees  of  intensity  of  the 
underlying  ionization.  For  each  such  value  of  critical 
frequency,  simulated  ionogram  generation  (synthesis)  was 
used  to  compute  the  values  of  minimum  group  paths,  again 
at  frequencies  of  6.5,  7.0,  and  7.5  MHz.  Denoting  these 
values  by  P'(6.5),  P'(7.0),  and  P'(7.5),  it  was  then 
possible  to  demonstrate  (as  shown  in  Figure  9)  the  signifi- 
cance of  underlying  ionization  upon  the  measurement  of 
minimum  group  path  through  the  use  of  an  RMS  error  defined 


as : 
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As  would  be  expected,  the  larger  values  of  critical  frequency 
in  the  underlying  layer  were  associated  with  larger  RMS 
errors.  Furthermore,  when  the  error  at  each  individual 
frequency  is  plotted  (as  shown  in  Figure  10)  it  can  be  seen 
that  the  effects  of  underlying  ionization,  at  each  of  the 
sounding  frequencies,  is  quite  similar  to  the  composite 
effect  shown  in  Figure  9. 

In  conclusion,  therefore,  our  example  demonstrates 
that  the  presence  of  an  underlying  layer  can,  indeed, 
affect  the  measurement  of  minimum  group  path,  and  that 
the  degree  to  which  the  measurement  is  affected  depends 
upon  the  "strength"  of  the  underlying  layer.  While  this 
conclusion  is  hardly  surprising,  it  is  significant  to 
note  that  even  when  realistic  estimates  were  used  to  model 
the  underlying  ionization,  the  measured  values  of  minimum 
group  path  were  substantially  altered  from  the  values  which 
would  be  measured  in  the  absence  of  such  ionization. 

A related  problem,  then,  is  to  determine  how  this 
alteration  in  the  group  path  will  affect  the  ionostate 
which  is  computed  by  the  inversion  program.  Again,  we 
have  proceeded  to  study  this  problem  through  the  use  of 
a particular  example  in  which  the  ionosphere  was  assumed 
to  consist  of  two  layers  — an  LQP  model  layer  and  a layer 
of  underlying  ionization.  The  LQP  model,  as  before,  was 
specified  by  the  ionostate  {6570,  6720,  5,  0,  0,  2},  while 
the  underlying  layer  was  set  to  have  a critical  frequency 
of  .75  MHz.  Actually,  the  electron  density  profile  in 


Figure  10.  Error  introduced  b 
sounding  frequency 
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this  case  is  given  by  Figure  8.  Our  next  step  was  to 
obtain  the  results  of  the  simulated  backscatter  ionogram 
generation.  Again,  these  results  were  available  from  some 
of  the  previous  runs  which  were  used  in  the  construction 
of  Figures  9 & 10,  and  are  summarized  below 

frequency  elevation  (degrees) min,  group  path  (km) 


33.89 

1081.4 

28.08 

1260.5 

The  next  step,  of  course,  was  to  use  the  values  of  minimum 
group  path,  given  in  the  table  above,  as  input  for  the 
inversion  program.  For  the  purpose  of  comparison,  it  is 
worthwhile  to  reconsider  Figure  9,  at  this  point.  We 
know  that  the  correct  ionostate  for  the  main  (LQP)  layer 
was  {6570,  6720,  5,  0,  0,  2}.  However,  when  the  underlying 
ionization  is  present,  we  can  see  that  this  ionostate  will 
result  in  a RMS  error  of  about  5 km  (or  a mean  squared 
error  of  25  km2).  This  error  then  provides  somewhat  of  a 
"tidemark"  in  the  sense  that  when  a lower  mean  squared 
error  is  produced  in  the  course  of  running  the  inversion 
program  for  this  particular  case,  it  would  seem  feasible 
that  the  program  is  attempting  to  construct  a composite 
(or  compromise)  ionosphere  model  to  simulate  the  effects 
of  both  the  main  and  underlying  layers  of  ionization. 
Indeed,  this  situation  was  found  to  occur.  The  following 
graph  illustrates  the  convergence  of  the  mean  squared  error 


Figure  11.  Two  component  (GFC&RB)  convergence  of  MSE  in 
presence  of  underlying  layer. 
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when  RB  and  GFC  are  specified  as  variable  ionostate  com- 
ponents. The  mean  squared  error  after  one  iteration  is 
approximately  7 km2,  as  opposed  to  the  25  km2  error  which 
would  be  obtained  for  the  apparently  correct  ionostate  of 
{6570,  6720,  5,  0,  0,  2}.  A second  iteration  was  also 
used,  but  failed  to  provide  a significant  change  in  either 
the  ionostate  or  the  mean  squared  error. 

2.4  QUADRATIC  GRADIENT 

In  this  section,  an  additional  ionostate  component 
is  introduced  in  order  to  allow  for  the  possibility  of  a 
"quadratic  gradient"  in  the  critical  frequency.  While  the 
resulting  model  may  still  be  regarded  as  locally  quasi- 
parabolic (LQP) , the  dependence  of  the  critical  frequency 
upon  the  geocentric  angle  subtended  by  the  ray  trajectory 
is  now  of  the  form: 

f c ( 0 ) = FC  + GFC*0  + G2FC*0  2 (20) 

with  GFC  and  G2FC  as  constants.  The  other  quasi-parabolic 
parameters  contain  a "linear  gradient",  i.e. 

rb(0)  = RB  + GRBx0 

r (0)  = RM  + GRMx0 

m 

Thus  if  we  assume,  as  was  previously  done,  that  the  over- 
head values  of  the  quasi-parabolic  parameters  are  known, 
the  present  section  is  intended  to  consider  the  simultaneous 
determination  of  four  remaining  ionostate  components  (i.e. 

GRB , GRM,  GFC,  and  G2FC) . In  order  to  test  this  capability 
in  the  program  we  proceeded,  as  in  previous  cases,  by 
selecting  ionostate  component  values  of: 

I a 
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RB  = 6570  km 
RM  = 6720  km 
FC  = 5 MHz 
GRB=  25  km/rad. 
GRM=  20  km/rad. 
GFC=  2 MHz/rad. 
G2FC=  16  MHz/rad.2 


and  simulating  the  measurement  of  a backscatter  ionogram 
in  an  LQP  model  ionosphere  at  sounding  frequencies  of 
6.5,  7.0,  7.5,  and  8.0  MHz.  The  results  of  this  simulation 
are  summarized  in  the  table  below: 

SOUNDING  FREQUENCY  (MHz.)  MINIMUM  GROUP  PATH  (km) 


6.5 

7.0 

7.5 

8.0 


1065.2 

1149.8 

1234.6 

1322.0 


These  values  of  minimum  group  path  were  then  used  as  part 

of  the  data  input  file  for  the  inversion  program.  In  addi- 

\ 

tion,  an  initial  ionostate  of 

RB  = 6570  km 

RM  = 6720  km 

FC  = 5 MHz 


GRB 


30  km/rad. 
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GRM  = 

28 

km/rad. 

GFC  = 

1. 

8 MHz/rad. 

G2FC  = 

13 

MHz/rad . 2 

was  pecified.  Allowing  the  inversion  program  to  run  on 
the  basis  of  this  data  resulted  in  the  generation  of  new 
estimates  for  the  values  of  the  variable  ionostate  com- 
ponents, as  shown  in  Figure  12.  Referring  to  this  partic- 
ular figure,  it  can  be  seen  that  the  final  set  of  variable 
ionostate  components  is  still  rather  far  away  from  the 
correct  value,  despite  the  fact  that  the  mean  squared 
error  is  reduced  to  . 4 km2 . There  are  two  possible  explan- 
ations for  this  paradoxical  result.  The  final  set  of  iono- 
state components  may  represent  an  additional  local  minimum 
in  the  mean  squared  error  hypersurface.  Alternatively, 
the  small  value  of  mean  squared  error  at  the  final  set  of 
ionostate  components  may  represent  the  offsetting  effects 
of  certain  combinations  of  variable  ionostate  components. 
More  specifically,  if  a backscatter  ionogram  leading  edge 
is  synthesized  in  both  the  correct  LQP  model  ionosphere 
and  the  final  LOP  model  ionosphere,  the  results  as  shown 
in  Figure  13  are  virtually  indistinguishable.  That  there 
is  any  deviation  at  all  becomes  apparent  when  the  difference 
in  minimum  group  paths  is  plotted  as  a function  of  frequency 
(Figure  14) . In  order  to  demonstrate  why  we  believe  that 
mutually  offsetting  effects  may  be  responsible  for  the  small 
value  of  me.  squared  error  at  the  final  set  of  ionostate 
components,  we  have  plotted  (in  Figure  15)  the  over-head 
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Figure  13. 


Backscatter  leading  edges  for  two  different 
ionostates . 


Fiqure  14.  Difference  in  backscatter  leading  edges  for 
two  different  ionostates. 
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Figure  15.  Critical  frequency  variation  for  two  different 
ionostates . 
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values  of  critical  frequency  as  a function  of  ground  range 
distance  from  the  sounder  location  (in  the  direction  of 
the  critical  frequency  gradient) . The  lower  curve  cor- 
responds to  the  final  set  of  ionostate  components  and  is 
a graph  of 

fc=5+ (2 . 28478) x (ground  range/6370) +(13.86538) x (ground  iange/6370)2 

(21) 

The  upper  curve  corresponds  to  the  correct  set  of  ionostate 
components  in  which  case 

fc  = 5+ (2) x (ground  range/6370)  + 16x (ground  range/6370) 2 

(22) 

The  initial  separation  of  the  two  curves  is  on  the  order 
of  .02  MHz  and  increases  only  to  about  .1  MHz  at  rather 
large  values  of  ground  range  (^2000  km) . In  a relative 
sense,  this  latter  divergence  represents  less  than  a 2% 
deviation  from  the  nominal  value,  and  thus  we  believe  that 
the  final  set  of  ionostate  components  was  strongly  influenced 
by  the  offsetting  effects  of  the  GFC  and  G2FC  components, 
or,  in  other  words,  we  believe  that  it  is  difficult  to 
resolve  the  critical  frequency  gradient  into  separate 
"linear"  and  "quadratic"  terms  when,  as  in  the  present 
case,  the  ground  ranges  are  small  and  when  the  "linear" 
and  "quadratic"  terms  contribute  in  a nearly  equal  manner 
to  the  overall  gradient  in  the  critical  frequency. 


We  did,  however,  consider  the  possibility  of  detecting 
a "quadratic"  component  in  the  absence  of  a "linear"  term. 
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In  this  case,  the  correct  ionostate  was  assumed  to  consist 
of 


RB  = 

6570  km. 

RI1  = 

6720  km . 

FC  = 

5 MHz 

GRB  = 

25  km/rad. 

GRM  = 

20  km/rad. 

GFC  = 

0 MHz/rad. 

G2FC= 

16  MHz/rad 

The  synthesis  program  was  then  run  in  order  to  compute  the 
minimum  group  paths  occuring  at  6.5  and  7.0  MHz.  These 
results  were  then  used  in  the  inversion  program  with  esti- 
mated initial  values  for  RM  and  G2F C.  The  convergence  of 
these  parameters  is  shown  in  Figure  16  for  two  sets  of 
initial  estimates.  It  can  be  seen  that  a rather  substan- 
tial reduction  in  the  mean  squared  error  has  occured. 
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3.  PROGRAM  DOCUMENTATION 

3.1  OVERVIEW  OF  PROGRAM  OPERATION 

The  purpose  of  this  section  is  to  provide  an  overall 
view  of  the  operation  of  the  backscatter  ionogram  inversion 
and  synthesis  programs.  Figure  17  demonstrates  how  the 
synthesis  portion  of  the  program  is  accomplished.  Since  the 
program  automatically  assumes  that  the  ionospheric  electron 
density  profile  follows  the  LQP  model,  only  the  values  of 
the  six  ionostate  components,  the  sounding  frequencies,  and 
a few  control  parameters  (which  will  be  discussed  later) 
need  to  be  entered.  The  program  then  computes  the  minimum 
group  paths  for  each  of  the  sounding  frequencies. 

The  inversion  operation,  as  shown  in  Figure  18  is  some- 
what more  complicated.  In  this  case,  the  input  consists 
of  an  estimated  ionostate,  a set  of  sounding  frequencies, 
the  corresponding  minimum  group  paths  which  were  observed 
at  each  sounding  frequency,  and  some  control  parameters. 

In  its  simplest  terms,  the  program  then  attempts  to  modify 
the  estimated  ionostate  in  a manner  which  results  in  the 
minimization  of  the  mean  squared  error.  It  should  be  noted, 
however,  that  while  Figures  17  and  18  illustrate  the  function 
of  the  program,  the  actual  details  of  the  program  flow  are 
considerably  more  complex.  It  is,  therefore,  our  intention 
to  provide  a more  detailed  discussion  at  this  time.  This 
discussion  will  be  limited  to  the  inversion  program  since 
the  synthesis  program  is  actually  contained  in  the  inversion 
program,  and  will  also  be  limited  to  the  major  subroutines. 
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Figure  17 


Synthesis  program  flow  chart 
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INPUT: 
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Figure  18.  Inversion  program  flow  chart. 
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3.2  FLOW  OF  MAJOR  SUBROUTINES 

The  most  fundamental  part  of  the  inversion  program  is 
the  ray  tracing  subroutine,  RAYT.  Referring  to  the  flow 
chart  provided  in  Figure  19,  the  first  step  upon  entering 
RAYT  is  to  increase  the  ray  counter  by  1.  The  quantities 
Y ( 1 ) through  Y(6)  then  correspond  to  the  terms  on  the  left 
hand  side  of  equations  1 through  6,  and  are  initialized  to- 
gether with  the  "trouble  indicator"  N(3).  Because  the 
initial  region  encountered  by  the  ray  is  free  space,  a 
simple  straight  line  ray  tracing  subroutine  (RAYS)  is  called 
to  trace  the  rays  from  the  source  location  (assumed  to  be 
on  the  surface  of  the  Earth)  to  the  base  of  the  ionosphere. 

At  this  point,  a loop  is  initiated  which  stores  the  present 
values  of  Y(l)  through  Y(6)  in  YY(1)  through  YY(6),  the 
latter  array  serving  as  a "place  keeper" . Upon  exiting 
this  loop  the  subroutine  consists,  in  essence,  of  a fourth 
order  Runge-Kutta  scheme  for  the  solution  of  the  ray  equa- 
tions. The  counter,  NSTP,  is  used  to  keep  track  of  the 
number  of  integration  steps  which  have  been  performed.  At 
this  point,  the  ray  tracing  program  consists  of  four  branches, 
labelled  M=1  through  M=4  and  corresponding  to  the  four  parts 
of  an  integration  step  in  the  4th  order  Runge-Kutta  solution. 
At  statement  200  a test  is  performed  to  guarantee  that  the 
ray  point  lies  within  the  spatial  domain  of  the  LQP  model 
i.e.  r^fr <r  . If  r<rb  the  ray  is  assumed  to  be  returning 
to  Earth  in  which  case  the  free  space  ray  tracing  subroutine 
RAYS  is  called.  If  r>rm,  it  is  assumed  that  the  ray  will 
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Figure  19  continued 
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Figure  19  continued. 
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never  return  to  Earth  and  therefore  the  trouble  indicator 
N(3)  is  set  to  1.  Most  of  the  time,  however,  the  ray  will 
be  in  the  spatial  domain  of  the  LQP  model  and  therefore  the 
right  hand  sides  of  the  ray  equations  1 through  6 will  be 
evaluated  and  stored  in  G(M,1)  through  G(M,6).  When  all 
four  branches  of  this  process  have  been  completed  (M=l 
through  M=4)  one  complete  integration  step  will  have  been 
completed  at  which  time  the  program  will  check  again  to 
see  if  it  should  continue.  Upon  returning  from  subroutine 
RAYT  the  program  will  either  have  a trouble  indicator  N(3) 
set  to  1 or  else  it  will  have  the  group  path  corresponding 
to  the  input  frequency  and  elevation  angle. 

This  information,  however,  is  not  sufficient  since 
we  are  interested  only  in  rays  which  have  the  minimum 
group  paths  for  a given  frequency.  For  this  reason, 
subroutine  RAYT  is  nested  inside  another  subroutine  (STAR) 
which  varies  the  initial  ray  elevation  angle  until  a minimum 
group  path  is  obtained. 

Subroutine  STAR  was  written  to  compute  the  minimum 
group  path  for  a ray  (with  a given  frequency)  traced 
through  an  ionosphere  with  given  parameters.  The  only 
parameter  which  may  be  varied  in  order  to  compute  this 
minimum  is  the  elevation  angle  of  the  ray. 

Referring  to  the  flow  chart,  an  initial  elevation 
angle  is  part  of  the  input  to  this  subroutine.  Subroutine 
RAYT  is  then  called  in  order  to  ensure  that  a ray  with  the 
initial  elevation  angle  will,  in  fact,  reflect  back  to  the 
surface  of  the  Earth.  If  this  situation  does  not  occur, 


yy 


Figure  20.  Flowchart  for  subroutine  STAR. 
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RAYT  will  set  the  flag  "N(3)=l"  and  thus  cause  RESTAR  to 
be  called. 

At  this  point,  assuming  that  RESTAR  has  not  been 
called,  the  initial  elevation  angle  is  saved,  the  ray 
counter  N(4)  is  initialized,  and  the  N(3)  flag  is  set 
to  zero.  The  program  then  constructs  a "star"  or  "cluster" 
which  consists  of  3 equi-spaced  elevation  angles  (denoted 
by  6(1),  6(2)  and  6(3))  (see  Fig.  21).  The  center  of  the 
star  is  the  initial  elevation  angle.  For  each  of  the  3 
elevation  angles  in  the  star,  a corresponding  group  path 
is  computed.  In  the  event  that  one  (or  more)  of  the  eleva- 
tion angles  will  not  produce  an  Earth-striking  ray,  the 
width  of  the  star  is  decreased  by  50%  and  the  process  is 
repeated. 

Assuming  that  all  3 rays  will  land  on  the  Earth's 
surface,  and  assuming  that  the  ray  counter  N ( 4 ) is  less 
than  N(5),  the  next  task  of  the  subroutine  is  to  determine 
if  the  value  of  group  path  at  the  center  of  the  star  is 
less  than  the  group  path  at  the  extremities.  When  this 
situation  occurs  (i.e.  G(1)>G(2)  and  G ( 3 ) >G ( 2 ))  the  minimum 
group  path  has  been  "localized".  Program  flow  is  then  sent 
to  statement  4 on  the  flow  chart.  This  action  initiates  a 
two  stage  "flatness  test".  The  first  stage  consists  of 
determining  if  either 

| G ( 3)  - G (2 ) | < S (83) 
or 

|G(1)  - 


G (2)  | < S ( 83) 
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is  satisfied.  If  neither  inequality  is  satisfied,  the 
star  size  is  decreased  by  50%  and  the  program  flow  is 
sent  back  to  statement  8.  However,  when  either  of  the 
inequalities  is  satisfied,  the  second  stage  flatness  test 
is  initiated.  In  beginning  this  process,  the  subroutine 
"realizes"  that  the  minimum  group  path  occurs  in  either 
the  interval  [8(1),  8(2)]  or  [8(2),  8(3)]  because  it  is 
assumed  that  the  group  delay  vs.  elevation  angle  curve  is 
concave.  It  is  also  true,  at  this  point  in  the  subroutine, 
that  the  difference  in  group  path  values  corresponding  to 
the  end  points  of  the  interval  is  less  than  S(83).  This, 
however,  is  not  always  sufficient  to  establish  the  minimum 
value  of  group  path,  as  can  be  seen  from  the  following 
example : 


I 
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Subroutine  steep  is  the  most  important  subroutine  in 
the  inversion  program.  As  can  be  seen  in  the  accompanying 
flow  diagrams  (Fig.  22),  this  subroutine  begins  by  computing 
a search  direction  which  corresponds  to  the  direction  of 
steepest  descent  along  the  mean  squared  error  hypersurface. 

NM  represents  the  maximum  tv’mber  of  allowable  search  direc- 
tions and  if  we  assume  that  the  current  number  of  search 
directions  is  less  than  this  limit,  the  subroutine  then 
computes  how  far  the  individual  ionostate  components  may 
be  varied  along  the  search  direction  before  the  search 
goes  outside  of  the  ionostate  domain.  As  an  example,  we 
may  consider  the  situation  shown  in  Fig.  23  in  which  the 

A 

unit  vector,  t,  is  assumed  to  represent  the  search  direction. 
The  two  variable  ionostate  components  denoted  here  by 
and  £2  are  assumed  to  be  in  the  domain  ^lMAX^  x 

^2MIN  ^2MAX^  * Now'  the  distance  from  to  the  line 

Climax  is'  of  course»  I ^imax“^10  I • Therefore,  the  dis- 
tance, s-^,  which  we  must  proceed  along  the  search  direction, 

A 

t,  in  order  for  to  assume  the  value  of 

I ^lMAX~^lol/  1^1  I * Similarly,  in  considering  £2,  we  find 

that  the  maximum  distance,  S2»  along  the  search  direction 

is  I | / 1 t9  | . Therefore,  in  this  case,  the  maximum 

2MAX  20  z 

A 

allowable  distance  along  the  search  direction,  t,  is  given 
by  s2  since  s2  is  less  than  s^.  These  methods  can  be 
readily  extended  to  the  simultaneous  consideration  of 
several  variable  ionostate  components,  as  shown  in  the 
accompanying  flow  chart.  Once  the  search  direction  and 
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Figure  22.  Flow  chart  for  subroutine  STEEP. 
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maximum  distance  have  been  established,  the  minimization 
procedure  be'comes  a one  dimensional  problem  in  which  our 
immediate  objective  is  to  find  an  ionostate  which  is  on 
the  search  line  and  results  in  minimum  mean  squared  error. 
This  one  dimensional  minimization  uses  the  same  scheme  as 
was  used  in  subroutine  STAR,  including  the  same  sort  of 
tests  which  were  used  for  termination  of  the  minimization. 

In  the  present  subroutine,  however,  the  fact  that  a minimum 
has  been  found  does  not  necessarily  mean  that  the  subroutine 
has  finished  its  task.  In  particular,  if  the  mean  squared 
error  is  still  larger  than  a value  specified  in  the  input 
file,  the  subroutine  will  go  back  and  compute  a new  search 
direction  corresponding  to  the  local  properties  of  the 
mean  squared  error  hypersurface  in  the  neighborhood  of 
the  present  ionostate. 


This  completes  the  discussion  of  the  major  subroutines. 
The  minor  subroutines  are  much  simpler  and  are  summarized 
in  the  following  table. 


Table  1:  MINOR  SUBROUTINES 


(SUBROUTINE 


CALLED  FROM 


CALLS  TO 


PURPOSE 


•SUPER 


SYNTH 


MAIN 


SUPER 


SETD,  STEEP, 
SYNTH,  ERCK 


supervisory 


STAR 


simulation  of 
backscatter 
sounding  data 


SETD 


SUPER 


Reads  data  file 


RAYS 


RAYT 


Free  space  ray 
tracing 


RESTAR 


STAR 


RAYT 


Restart  minimum 
group  path  search 


ERCK 


SUPER, STEEP, 
GRADl 


STAR 


Compute  MSE 


GRADl 


STEEP 


ERCK 


Compute  gradient 
of  MSE 
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3.3  PROGRAM  LISTING 

The  following  pages  provide  a listing  of  the  present 
version  of  the  program.  The  file,  DF,  is  the  data  file. 
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SUBROUTINE  RAVI  74/74  OPT=1  TRACE 


FTN  4.6+439 


1 


5 


10 


15 


2 J 


25 


30 


35 


40 


45 


50 


55 


C 

C 

C 

C 


C 

C 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


c 

c 

c 


MIS  PAGE  IS  BEST  QUALITY  PRACTICABLE. 
*SOM  COPY  FURBISHED  TO  DDC 

SUBROUTINE  RAYT 
COMMON  S,Y,W 

Dj.MENSj.ON  S(200)  ,Y(6)  ,N  (200) 

DIMENSION  YY  (6)  ,G  (4,6)  ,SRX  (6) 

RL  = S (108) 

N (4)  =N  (4)  ♦ 1 

THE  PURPOSE  OF  THIS  SUBROUTINE  IS  TO  TRACE 
A RAY  OF  GIVEN  FREQUENCY  o ELEVATION  ANGLE. 

WHEN  RAY  IS  IN  FREE  SPACt.  SUBROUTINE  RAYS 
IS  CALLED. 

Y ( 1 ) = 6 3 7 0 . 

Y (2)  =0.  0 

Y 13)=S1N  ( S ( 1 ) ) 

1 (4)  =COS  (S  ( 1 ) ) 

Y (5)=0. 

Y ( 6 ) = 0 . 

N (3)  =0 

INITIALIZATION  COMPLETE,  RAY  IS  PRESUMABLY  IN 
FREE  SPACE.  THEREFORE,  RAYS  WILL  BE  USED 

CALL  RAYS 
DO  20  1=1, e 
20  YY  (I) = Y (I) 

NS I P=0 

THE  FOLLOWING  IS  A FOURTH  ORDER  RUNGE 
KUTTA  INTEGRATION 

1 1 M = 1 

GO  TO  200 

1 DO  2 1*1,6 

2 Y (I)  = YY  (I)  + (3  (20)/2.)  (!J,I) 

1ST  RUNGE  KUTTA  STEP  COMPLETED. 

M=  2 

GO  TO  200 

3 DC  4 1-1,6 

4 Y (I)=YY  (l)  + (S(20)/2.)  *G  (.1 ,1) 

2ND  RUNGE  KUTTA  STEP  COMPLETtD 
M = 3 

GO  TO  200 

5 DC  6 1=1,6 

6 Y (I)=YY  (i)+S(2Q)*G(M,I) 

3UL  RUNGE  KUTTA  STEP  COMPLETED 
ft*  4 


THIS  PACE  IS  BEST  QUALITY  PRACTICAL 
FROM  OOPY  FURNISHED  TO  DDC 
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SUBROUTINE  R A Y 1 *74/74  OPI=1  TRACE 


FT N 4.6*439 


6 0 


7 J 


75 


80 


as 


9 0 


95 


100 


1)5 


110 


GO  TO  200 

7 DC  a 1=1,6 
SRX  (I)  =YY  (I) 

8 YY(l)=iY  (I)  ♦ (S(2  0)/b.)*  (O  (1,1)  *2.*G(2,I)  *2.  *G  (3 , 1)  ♦ G (4 , 1)  ) 

C 

C FINAL  LJNGL  KUTTA  STEP  COMPLETED 

C 

N STP=N  SI P* 1 
C 

C CHECK  FOR  NUMBER  OF  STEPS 

C 

IF  (NSTP.GT.  N (1)  ) GO  TO  101 
C 

C CHECK  FOu  FELL  SPACE 

tc  13= S (21)  +YY  (2)  *S  (24) 

If  (YY  (l).LT.S  (10HJ)  GO  TO  27 
GO  TO  11 
1 ) 1 CONTINUE 
PHI H T 102 

102  FORMAT  (1u  , 9HNSTP>N (1 ) ) 

STOP 

C 

C THIS  COMPLETES  THE  RUNG  L KUTTA 

C PART  OF  SUBROUTINE.  THE  FOLLOWING  PROCEDURE 

C IS  USEL  TO  EVALUATE  RIGHT  HAND  SIDE 

C OF  RAY  EQUATIONS. 

C 

200  RB=S  (21)  +Y  (2)  *S  (24) 

IF  (Y  (1)  .GE.RL)  GO  TO  215 
27  DC  21t>  K V=  1 , 6 
C 

C HACK  UP  RAY  TRAJECTORY  BY  ONE 

C STEP  TO  INITIALIZE  FREE  SPACE 

C RAY  TRACING  SUBROUTINE. 

C 

216  Y (KV)=SHX (KV) 

CALL  RAYS 
RETURN 

215  CONTINUE 

RM=S  (22)  *Y  (2)  *S  (25) 

IF  (Y  (1)  .LE.Rfl)  GO  TO  218 
C 

C CHECK  FOR  PENETRATING  RAY 

C 

N (3)  = 1 
S (27)  * 1000. 

RETURN 

218  FC=S  (23)  ♦ Y (2)  *S  (26) 

C 

C THE  FOLLOWING  STATEMENTS 

C ARE  USEL  TO  EVALUATE  "FORCING  FUNCTION" 

C TERMS  ON  THE  RIGHT  SIDE  OF  THE 

C RAY  EQUATION. 

C 

F=S  (28)  /PC 
YM=RM-RB 

T ER M 1=  1 . - (1 ./  (P**2)  ) 
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THIS  PASS  IS  BEST  QUALITY  PRACTICABLE 
FROM  COPY  FURNISHED  TO  DDC 


I 

| SUBROUTINE  RAY!  74/74  0?T=1  TRACE 


FT N 4.6*439 


115 


12  0 


125 


1 JO 


135 


140 


Ttl>»2=  ( (RM-Y  (1;  ) / (F*YM)  ) * (RB/Y  (1)  ) 

TERM2=TERM2**2 
Zi,lA=T  ERM  1 ♦ ! ERM2 
IF  (ZF.TA.GT.  1.00)  00  TO  76 

PY  2=Ziil  A 
PY=SQRT  (ZLT  A) 

£ I A=  (KF.-Y  (1)  ) / (F*YM) 

GAMMA=  (hii/Y  ( 1)  ) **2 

JPDR=  (-2.*ETh*GAMMA/(F*  YM)  > -2.  * (ETA**2/Y  ( 1)  ) * GAMMA 
DPDR=DPDR/(2.*PY) 

El A2  = ETA  * * 2 

DPD FC=  (1  ./PY)  * ( (-FC/(S  (28) **2) ) *ET A2*GAMM A/FC) 

DPDRM= (l./rY) * (GAMMA* ETA*  ( (GAMM A/  (HM-Y  ( 1) ) ) -GAMMA/YM) ) 
DPDRB=  (1./PY) *( (G AMMA*ET  A 2/YH)  ♦ (G AM M A*ET A2/RB ) ) 
DPEA=DPDRB*S  (24)  *DPDRM*S  (25)  +DPDFC*  (S  (26)  ) 

GO  TO  77 

76  PY=1.0 
PY2=  1. 0 
DPD3=orCA=0.0 

77  G (K,  1)  = Y (3)  /PY2 
G(M,2)=Y(4)/(Y(1)*PY2) 

G (M,3)  = (DPL'R/P Y)  ♦ (Y  (4 ) **2) / (Y  ( 1 ) *PY2) 

G (M  , 4)  = ( LPDA/  (Y  (1)*PY))  - ( (Y  (4 ) * Y (3)  ) / (Y  ( 1 ) * PY  2j) 

G (M,5)  =1./PY2 

G (M  , 6)  =1.0 

IF  (M.  Ew-  1)  i>0  TO  1 

IF  (1.  LQ.  2)  GO  TO  3 

lF(M.EU.-3)  GO  TO  5 

GC  TO  7 

END 


. ™ t c best  QUALITY  mCIICABI* 
JS  COPY  PUNISHED  TO  
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SJ3K0TITINE  SITE  74/74  ORT=1  TRACE  FTN  4.b«-439 


1 C 

C 

c 

c 

5 SUBROUTINE  SETD 

COMMON  S,Y,N 

DIMENSION  S (200)  ,Y  (6)  ,N  (2C0) 
C 

C NULL  OU'I  S 6 N ARRAYS 

1 J C 

LC  1 1=1, 2C0 

s (i)  =u.o 

1 N (!)  =0 

r 

\rj  i_  Not.  RE/iD  S S N ARRAYS 

C 

20  CONTINUE 

REAL'  (10,3)  K,Z 

3 FORMAT  (13,  IX,  E20.  7) 

?.J  IF(K.Lt.O)  GU  TO  4 

S (K) =Z 
GO  TO  20 

4 CONTINUE 

Fi-IAU  111,6)  K , M 

26  ( FORMAT  (13,15) 

IF  (K.  LF  . 3)  JO  TO  7 
N (K)  =M 
GO  TO  4 
7 RETURN 

3 J END 


SUSP.  OUT  INS  RAYS 


1 


5 


1 J 


15 


2 J 


2 5 


10 


j 5 


'4  0 


4 5 


50 


55 


C 

C 

c 

c 

SUBROUTINE  RAYS 
COMMON  S,Y,N 

DIMENSION  S (200)  ,Y  (b)  ,N  (2C0) 
HFP1=  3.  14150/2. 

P H 1=  S ( 1 ) 

Rs=S  (31) 

J 33  FORMAT  (111  , J (E2 0. 7,  1 X) ) 

C 

C THIS  SUBROUTINE  IS  CALLED 

C ONLY  WHEN  TilE  RAY  IS  IN 

C FREE  SI  ACE.  THIS  SUBROUTINE 

C COM PUT I S A STRAIGHT  LINE 

C TRAJECTORY. 

C 

if  (Y  (J)  . Lt.  9.  0)  GC  TO  4 
C 

c:  This  DETERMINES  IF  A AY  IS 

L GOING  Up  UR  DOWN 
C 

GC  TO  J 

3 RX=S  ( 1 Do) 

C 

C u AY  IS  GOING  UP,  BUT 

C HE  IS  CONSTANT 

C 

L 1-RF* Y (4)  /i(X 
TB1  = ASIN  (Z  1) 

THET  A - li  FPI-  Phi-  Til  1 
ST  EST=  IFl/ZI) *51 N (THETA) 

C 

C Inis  COMPUTES  THE  JPG01NG 

C PART  OF  THE  SUBROUTINE 

C 

1 i (1) =RX 

Y (2)  - ill  ETA 

1 h 1 = HE  PI-Tli  1 

Y (3)  =S1N  (TH  1) 

Y (4)  =COS  (Th  1 ) 

Y (5)  = Y (6)  -STE5T 
RETURN 

4 CONTINUE 
C * 

C THIS  PART  IS  FOR  DOWN  GOING 

C RAYS. 

C 

PHOOE Y- ACOF (hE/Y (1)  ) 

E2TA=ALS  (Y  (4)  ) 

ABPHI  = AC05  (SETA) 

1?  (ABPiii..LT.P*iOOEY)  GO  TO  5 
C 

C THIS  CLICK  IS  FOR  RAYS  WHICH 

C DC  NOT  CESCi.NL 

C STEEPLY  ENOUGH  10  TOUCH  THE 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
FROM  COPY  FURNISHED  TO  DD,Q  — 


SUBROUTINE  PAYS  74/74  0?T=1  TRACt  FT  N 4. 6*439 

C SURFACE  OF  lllx.  lAETH. 

C 

bO  A F L = l.  F P 1-  A B P ti  1 

AH  = ASIi!)  (Y  (1)  *3IN  (Ahji)/KE) 

AY  1 = 2.  J*iiFPi-AY  1 
ASTEST= (2.*nFPl) -AY 1-ARL 
SiLST=SIN  (ASTL5T)  *Y  (1  )/SIN  (A Y 1 ) 
b 5 Y (5)  = Y (5)  ♦SliST 

Y (6)  = Y (6)  + S'iE5T 

Y (2)  = Y 1 2)  ♦ASIEST 

Y ( 1 ) = P.  L 
S (27)  =Y  (5) 

70  RETURN 

S N (3)  = 1 

Y (5)  = Y (b)  = 1.  3^*03 

S(27)=Y(5) 

RETURN 

75  END 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
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SJBRCUTINE  ERCK  74/74  OPT=1  TRACE 


FT M 4.6  + 439 


1 


5 


10 


15 


2 J 


2 5 


3 3 


J 5 


40 


C 

c 

c 

c 

3UDROUIII1E  EiiCK 
COMMON  S , Y , N 

LI  ME  NS  ION  S (200)  , Y(6)  ,N  (200) 

PRINT  333 

333  FORMAT  ( 1 H , 4 HZRC K) 

C 

C This  SUBROUTINE  COMPUTE^  TiiE 

C At.AU  SQUARED  ERROR. 

C 

NF=N  (2) 

S (33)  =0 . 

LO  1 1=1,  MF 
C 

C s’flEQU£NCY  LOOP 

c 

I 1 = 1 + 7 
12=1+1 
13=1+13 
S (1)=S  (13) 

S (28) =S (12) 

C 

C SIT  FStjULNCY  8 INITIAL  ELEVATION  ANOLE, 

C CALL  STAR  TO  FIND  MINIMUM  OROUP  PATH. 

C 

CmLL  star 
5(I3)=S(1) 

AN3=1  J0.0*S(1)/3.  14  159 

1 S (33)  =3  (3  3)+  (((S(2  7)-5(I1))**2)/N?) 

IF  (N  (21)  .Sg.O)  30  TO  2 

PRINT  3',  (S(L)  , L=  2 1,  26) 

3 FORMAT ( 13  , E20. 7) 

PRINT  4 ,S  (33) 

4 FORMAT  (111  ,4ntlSE=,  U,  E20.  7) 

STOP 

2 RETURN 
END 


THIS  PACE  IS  BEST  QUALt TY PRACTICABXJ 

PROM  COPY  FURNISHED  TO  DDC  ' 70 


SUBROUTINE  STAR  74/74  OPT=1  TRACE  FT  N 4.6*439 


1 


5 


10 


13 


20 


2 3 


3 3 


35 


4 3 


45 


50 


C 

C 

C 

C 

SUBROUTINE  STAR 
COMMON  S,Y,N 

DIMENSION  S (200)  ,Y  (6)  fN  (200) 

C 

C ThlS  SUBROUTINE  FINDS  ThE  ELEVATION 

C ANGLE  WHICH  RESULTS  IN  MINIMUM 

C GROUP  PATH  AT  A GIVEN  FREQUENCY. 

C IT  ALSO  EVALUATES  THE  MINIMUM  GROUP 

C PATH. 

C 

DIMENSION  LhTA  (3)  ,U  (3) 

C 

C CHECK  ACCEPTABILITY  OF  INITIAL 

C ELEVATION  ANGLE. 

C 

CALL  RAYT 

IF  (N  (3)  . EQ.  1)  CALL  REST AR 
C 

C INITIALIZE  ENTRY  VALUE  OF  S(1) 

C 

S F = SL=S  (1) 

N ( 4)  = 1 
8 N (3)  =0 
C 

C SET  UP  STAR  6 EVALUATE 

C GROUP  PATHS 

C 

11  CONTINUE 
DC  1 1=1,3 

BETA  (I)  =SL*  (i-1j  *S  (84) 

S ( 1)  =BETA  (I) 

CALL  RAYT 

IF  (N  (3)  .SQ.0)  GO  TO  12 
PRINT  13 

13  FORMAT  (1H  , dOX,  1 3H*ST AR  S HR  INK  1 ) 

S ( 84 ) = S (o4)/2.  3 
GO  TO  11 

12  CONTINUE 

IF  (N  (4)  .GT.  N (b)  ) GO  TO  7 
1 G (I)  =S  (27) 

C 

C NOW  CHECK  TOR  MINIMUM  GROUP 

C PATH  AT  CENTER  OF  CLUSTER 

C 

J IF  (G  (1 ) .GT.G  (2)  ) GO  TO  2 
C 

C SHIFT  STAR  TO  SMALLEh  lLEVATICN 

C ANGLES. 

C 

G (3)=G  (2) 

G ( 2)  =G  (1) 

BET  A (3)  = SET  A (2) 

BhT  A (2)  = B ETA  (1) 
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S (1)=biTA  (1)  =BETA  (2)-S  (84)  FT N 4.6*439 

N 43) =0 
CALL  HA  YT 

IF  (N  (4)  .GT.N  (5)  ) GO  TO  7 
IF  (N  (3)  .EQ.  ))  GO  TO  14 
PRINT  15 

15  FORMAT  (111  ,6UX,  1 3li*STAR  SKR1NK2) 

S (84)  = S (d4)/2.0 

SL=  BETA (2) 

GO  TO  11 
14  CONTINUE 
G (1 ) =S  (27) 

GC  TO  3 

2 IF  (G  (3)  .GT.G  (2)  ) GO  TO  4 
C 

c shift  star  to  larger 

C ELEVATION  ANGLES. 

C 

G(1)=G(2) 

G (2)  = G ( 3) 

3ETA  (1)  =UETA  (2) 

BETA  (2)  = BETA  (3) 

S (1)=bc.TA  (3)  = BET  A (2)  *S  (84) 

N (3)  =0 
CALL  t:  A YT 

IF  (N  (4)  . Gl.  N (5)  ) GO  TO  7 
IF  (N  (3)  . cQ.O)  GO  TO  16 
PRINT  17 

17  FORMAT  (Hi  , BOX,  13b*STAR  SHRINKS) 

S (84)  =S  (84J/2.0 
SL=BET  A (2) 

GO  TO  11 

16  CONTINUE 

G (3)  =S  (27) 

GC  TO  2 
4 S L—  BE  1 A (2) 

C 

C THE  LEAST  GROUP  PATH  NOW 

C OCCURS  AT  HIE  CENTER  OF 

C TrlE  CLUSTER.  NO*  THIS 

C MINIMUM  VALUE  WILL  BE  TESTED. 

C 

T U=  ABS  (G  (3)  -G  (2)  ) 

TL=  ABS  (G  (1)  -G  (2)  ) 

IF  (TL.GT.S  (83))  GO  TO  5 
C 

C NOW  TEST  MIDPOINT  BETWEEN 

C B ET  A ( 1 ) 5 BETA  (2) 

C 

S (1)  = (uETA  (1)  *BST  A (2)  )/2.C 
N (3)  =0 
CALL  R A YT 

IF  (N  (4)  .GT.  N (5)  ) GO  TO  7 
IF  (A  (3 ) . Eg.  1 ) CALL  RESTAR 
E=  ADS  (S  (27)  -G  (2)  ) 

IF  (E. LE.S  (83) ) RETURN 
GC  TO  6 
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5 IF  (Id.  31.  S (83)  ) GO  TO  6 

NOW  TEST  MiUPOI NT  BETWEEN 
BETA  (2 ) 8 Ei.T  A ( 3 ) 

S (1)  = (BETA  (3)  + BET  A (2)  )/2. 

N (3) =0 
CALL  R A 71- 

IP  (N  (4)  .GT.  N (5)  ) GO  TO  7 
IF  IN  (3)  .EU.  1)  CALI  RESTAR 
E=  ABS (S (27) -G  (2) ) 

IF (E.LE.S  (b3) ) RETURN 

6 S (€4)=S  (84) /2.0 
SL=BETA (2) 

GO  TO  8 

7 CONTINUE 
PRINT  333 

333  FORMAT  (111  ,3UB,E) 

TOO  MANY  CALLS  TO  kAYT. 

TERMINATE  PROGRAM 

PRINT  21 

21  FORMAT  (Hi  , 8 OX, 1 6U* LAST  1CNOSTATE : ) 
PRINT  22,  (S(L)  , L = 21  ,2 6) 

22  format  ( i u ,aox,2n**,El5. 7) 

PRINT  23 

23  FORMAT  (1H  , 80 X , 1 8H* *F REQU E NCY  (MHZ):) 
PRINT  22, S (2d) 

PRINT  24 

24  FORMAT  (Iri  , 80X,  1 0U**CRASH00  1) 

STOP 

ENH 
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subroutine  iustau 

COMMON  5,Y,N 

DIMENSION  S (200;  , Y (i>)  ,N  (200) 
PRINT  1 

1 FORMAT  ( 1 ii  , 1 7;i*  RESTART  REQUIRED) 
K E=  0 
KMI N=0 

nl  PI=3  • 141*30/2.  J 
S 1KAX=  tiF  Tl~  1 . OE-06 
S 1MIN=  1.  DE-JO 
G M 1 N = 1 . 02  + Ub 

INITIALIZATION  FOR  FIRST  SEARCH 
13  COMPLETE. 

3  D5= (51 KAK-S1 MlN)  /30.  0 
DO  2 1=1,29 
S (1)=i*CS 

call  rayt 

WAS  KAY  SUCCESSFULLY  TRACED? 


IE  (N  (3)  . L\i.  1)  GO  TO  2 

IS  GROUP  PATH  > GMIN? 

IF  (S  (27)  .GT.GMIN)  GO  TO  2 
GMIN=S  (27) 

KMIN=1 
S 1 STAR=  S (1) 

2 CONTINUE 

Ii  AS  A ROUGH  MINIMUM  BEL  N FOUND?? 

IF  (KMIN.  F.Q.0)  GO  TO  5 

IS  SECOND  SEARCH  COMPLETE? 

IF  (KB. EQ.  1)  GO  TO  4 
S 1MAX=51STAR  + D5-1 . JE-06 
S1MIN=51STAR-DS+  1. J2- J6 
KB=  1 
GC  TO  3 

4 S ( 1)  = S 1 STAR 
RETURN 

5 CONTINUE 
PRINT  21 

21  FORMAT  (IB  , 80X,  17H**LAST  IONOST ATE: ) 
PRINT  22,  (S  (L)  ,L=21  ,26) 

22  FORMAT  ( 1 11  ,80X,2H**,E15.7) 

PRINT  23 

23  FOR  MAT  (III  , 30X,  1 8H**PBEyU  ENC  Y (MUZ):) 
PRINT  22,5  (28) 

PRINT  24 
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24  FOtwTAT  ( 1 il  , 8 OX,  1 0rt**CRASilC02) 
STOP 
END 
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C 

C 

c 

c 

SUBROUTINE  GRAD  1 
COMMON  S,  Y,  N 

DIMENSION  S (200)  , Y (6)  , N (2  00) 

Cl  MENS  ION  P2  (200) 

C 

C THIS  SUBROUTINE  COMPUTES 

C THE  GRADIENT  OF  THE 

C BEAN  SQUARED  ERROR. 

C FIRST,  THE  INITIAL  MEAN  SQUARED  ERROR 

C IS  STORED  IN  A SEPARATE  LOCATION. 

C 

CALL  ERCK 
S (66)  =S  (33) 

PRINT  366,  (S  (LL)  ,LL=14,  16) 

386  FORMAT  (1H  ,3  (E20. 7, IX ) ) 

EPS  1 = 1 0 . 0*S  (94) 

LPS2= 1 000. 0*S  (94) 

C 

C SET  GRADIENT  INDEX,  IG. 

C 

1G=4 1 
C 

C ENTER  SEARCH  LOOP  FOR 

C VARIABLE  ICNOST ATE  COMPONENT. 

C 

NTIST=G 
DO  3 1=1,6 
i N=  I * 5 
ID=I+b6 
IS=I*20 

IF  (N  (IN)  . EQ.  0)  GC  TO  3 
9 DS=  S (ID) 

IF  (DS.  EQ.  0)  GO  TO  4 
C 

C A VARIABLE  IONOSTATE  COMPONENT 

C HAS  BEEN  FOUND.  THE  PARTIAL 

C DERIVATIVE 'OF  MSE  WITH  RESPECT 

C TO  THIS  COMPONENT  WILL  NOW  BE 

C CALCULATED. 

C 

SCLD=3  (IS) 

1 3 S (IS)=S  (IS)  ♦ DS 
CALL  ERCK 
E 1 = S (33) 

S (IS)  = S (IS)  - (2.  0*DS) 

CALL  ERCK 
E2=S  (33) 

L AB= ABS  (E1-E2) 

C 

C CHECK  THE  DIFFERENCE  IN  MSE. 

C IP  IT  IS  TOO  SMALL,  IT  IS 

C VULNERABLE  TO  NUMERICAL  ERROR. 

c if  it  is  too  large,  it  is  not 
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C 
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C 

C 

C 


3 
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d5 


9 6 

L 
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28 

26 

27 

29 


C 

C 

C 


21 

22 

23 

<1 


166 


165 


25 
1C 
12 
14 
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REPRESENTATIVE  OF  LOCAL  PEOPLPTi.ES. 

IF  (LAB.  LI.  EPS  1)  GO  TO  16 
IF  (FAB.GT.EPS2)  GO  TO  11 

NOW  RESET  S(1S),  PACK 
Til E GRADIENT  VECTOR 
AND  CONTINUE  THE  SEARCH 
LCCP. 

S (IS)=30LD 

S (IG)  = (L1-£2)/(2.  6* DS ) 

P2  (I G)  = (L1*E2-2.0*S  (86)  )/  (4.*DS*DS) 

CHECK  FOR  LOCAL  CONCAVITY 

IF  (P2  (1G)  .LE.0.0)  NTEST=1 

IG=IG+  1 

CONTINUE 

NN=  4 1 + N ( 1 2) - 1 

IF (NTLST. EQ- 0)  GO  TO  26 

PRINT  28 

FORMAT  (lii  ,861,  7U*C0N  VEX) 

RETURN 

DO  27  K=  4 1 , ti  N 
3(K)=S(K)/P2(K) 

PRINT  29 

FORMAT  11H  , 8 0 X , 8 H *C  0 N C A V E ) 

RETURN 

FATAL  ERRORS 

FORMAT  4 1 H , 80X  , 1 7 ii*  *L  AST  IONOSTATE:  ) 
FORMAT  (1H  ,60X,2H**,E15.7) 

FORMAT ( IE  , 80X, 1 7H*  *MLAN  SQ.  ERROR: ) 
CO NT IN Ur 
PRINT  21 

PRINT  22,  (S  (L)  ,L=21  ,26) 

PRINT  23 
PRINT  22, S (66) 

PRINT  25 

FORMAT  (1H  , 8 OX , 1 OH* *C RASH 00 4) 

STOP 

DS=2. 0*DS 
PRINT  12, DS 

FORMAT  ( 1 h , 80X, 4H *DS= , L20 . 7) 

PRINT  14,1 

FO  RMAT  ( 1H  , 86X, 9H*3TATE  *-,13) 

GO  TO  13 
DS=DS/2.0 
PRINT  12, DS 
PRINT  14,1 
GO  TO  13 
END 
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1 C 

c 

c 

c 

5 SUBROUTINE  SUPER 

common  s,y,n 

DIMENSION  S (200)  , Y (6)  , N (2C0) 
DIMENSION  X (6)  ,G  (6)  ,H  (40) 

C 

10  C THIS  IS  A SUPERVISORY  SUBROUTINE 

C 

CALL  SETD 

IF  (N  (16)  .EQ.  1)  GO  TO  8 

IF  (N  (21)  . EQ.  1)  GO  TO  12 

15  CALL  STEEP 

DO  2 1=21,26 
PRINT  3,  S(I) 

3 FORMAT  '111  , E20. 7) 

2 CONTINUE 

20  PRINT  1 00, S (33) 

100  FORMAT (1H  , 4HMS E= , 1 X, E2 0. 7) 

STOP 

8 CALL  SYNTH 
STOP 

25  12  CALL  ERCK 

PRINT  1 0 1 ,S  (33) 

101  FORMAT  (Iri  ,411NSE=,E20.7) 
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1 C 

C 

c 

c 

5 SUBROUTINE  SYNTH 

COMMON  S,Y,N 

DIMENSION  S (200)  , Y (6)  , N (200) 

C 

C THE  PURPOSE  OF  THIS  SUBROUTINE 

13  C IS  TO  SYNTHESIZE  A BACKSCATTER 

C IONOGhAM. 

C 

NB=N  (2) 

DO  1 1=1,  ND 
15  S (1)  =S  (1*13) 

S (28)  =S  (1*1) 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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C 

c 

c 

c 


1 

24 

50 

51 


IOC 

101 

102 


SUBROUTINE  STEEP 

THIS  IS  A SPECIAL  STEEPEST 
CESCENT  SEARCH  SUBROUTINE  TO 
DETERMINE  OPTIMUM  IONOSTAT*. 

SEARCH  IS  CONSTRAINED  TO  DOMAIN 
DEFINED  BY  THE  INTERVAL  (ZL  (I)  , ZU  (I)  ) 

FOR  1=1,6  HH1CH  MUST  CONTAIN  THE 
ORIGINAL  ESTIMATED  I0N0S7 ATE. 

SEARCH  HILL  BE  TERMINATED  WHEN  ONE 
OF  THE  FOLLOWING  CONDITIONS  IS  FOUND 
TC  OCCUR: 

l.i  CALLS  TO  ERCK  • GT.  N (14) 

2. HSE.LT.S  (93) 

3. #  SEARCH  D1RECTICNS=N  (20) 

CCHMON  S, Y, N 

DIMENSION  S(2Q0)  , Y (6)  ,N  (200) 

DIMENSION  X (6)  , ZL  (6)  , ZU  (6)  , T (6)  ,B  (6)  ,H  (6) 
KOUNT=  0 
NM=N  (20) 

NE=N  (12) 

N (20)  = 0 
K=  1 

DETERMINE  SEARCH  DOMAIN 
STORE  10N0STA1E  IN  X 

DO  1 1=1,6 
IN=I*5 

IF  (N  (IN)  .NE.  1)  GC  TO  1 
ZL  (K)  = S (1*94) 

ZU  (K)  =S  (1*100) 

X (K)=S  (1*20) 

K=K*  1 
CONTINUE 
CONTINUE 
CALL  GRAD1 
RCUNT=KCUNT* 1 
GNORN=0.0 
DO  2 1=1,  NB 

GNORH=GNORM* (S(I*4J)  **2) 

GNORH=SCBT  (GNORM) 

DO  4 1=1,  NH 
T (I)  =~S  (1*40)  /GNORM 
PRINT  3 

FORMAT (1H  ,80X, 1JH* STATUS  OK) 

PRINT  100,N (20) 

FORMAT  (1H  ,bOX, 9 HIT ER AT ION, 13) 

PRINT  101 

FORMAT  (1H  ,80X,  18HPRJJSENT  IONOSTATE: ) 
PRINT  102,  (S(L)  ,L=2 1,26) 

FORMAT (1H  ,80X,S20.7) 
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c 

c 
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PRINT  103 

FORMAT  (1H  ,8  OX, 17  dSEARCH  DIRECTION : ) 
PRINT  102,  (T  (1) ,1=1,6) 

PRINT  104, S (86) 

FORMAT  1 1H  ,80X,4HMSE=,E20.7) 

IF  (N  (20)  .LI.  NH)  GO  TO  105 
PRINT  106 

FORM  AT  ( 1 H ,80X, 10H**CRASHC05) 

STOP 

ft  (20)  =N  (20)  *1 

NOW  DETERMINE  MAXIMUM  POSSIBLE 
DISTANCE  ALONG  THIS  DIRECTION. 

D=  (ZL  (1)-X  (1)  )/T(1) 

E—  (ZO(I)-X(I)  )/T  ( 1) 

ED=D*E 

IF  (ED.  GE.  0.)  GO  TO  5 
IF  (E.GT'«  0.  ) GO  TO  6 
E=D 

DM AX=E 

NOW  PERFORM  SAME  COMPUTATION  FOR 
REMAINING  DIMENSIONS. 

•DO  7 1 = 2, NB 

D=  (ZL(i)-X(l)  )/T(I) 

E=(ZU(i)-X(I))/T(I) 

EL=E*D 

IF  (ED.GE.0.)  GO  TO  5 
IF  (F..GT.0.)  GO  TO  8 
E=  D 

DMAX= AMIN  1 (DMAX, 2) 

CONTINUE 

INITIALIZE  STAR  SIZE 
AND  BEGIN  SEARCH 

DELT=DM AX/3. 0 
H (1)  =0. 

DX0=H (2) = D E L T 
11  (3)  =2.  0*DX0 
DO  9 J=  1 , 3 
K=  1 

DC  10  1=1, b 
1 N=I*5 

IF  (N  (IN)  . NE.  1)  GO  TO  10 
IS= 1*20 

S (IS)=X  (K)  *H  (J)  *T  (K) 

K=  K*  1 
CONTINUE 
CALL  ERCK 
KOUNT=XCUNT* 1 
B(J)=S(33) 

IF  (6(1j  . GT.  b (2)  ) GO  TO  11 
D (3)  =B (2) 

O (2)  =B  ( 1) 
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113 


123 


123 


113 


1 35 


143 


145 


133 
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1b  >3 


165 


170 


H (3)  =11  (2) 

DX0  = h (2)  =H  (1) 

IF  (TX  3. IT . 0 . ) GO  TO  5 
IF  (OX3.GT.LKAX)  GC  TO  5 
11  ( 1 ) = uXO-UF  LT 
K-  1 

DO  12  1=1, t 
I N=I*5 

IF  (N  (IN)  . NL.  1)  GO  TO  12 
10=1*20 

S ( 1 S)  = X (K)  * H (1)  *T  (K) 

K=K*  1 

12  GO NT 1 MU 1 
CALL  LRCK 
KOUNT=KCU!ll*1 
D(1)=S(j3) 

GC  TO  13 

11  IF  (D  (3)  .GT.B  (2)  ) GO  TO  14 
3(1)=B  (2) 

B(2)=B(3) 

a (i)  =ii  i2) 

bX3  = H (2) =H  (3) 

IF  (DXO. LT.O. 0)  GO  TO  5 
IF  (DXO.GT.LMAX)  GO  TO  5 
H (3)  =DXU*DILT 
K=  1 

DC  15  1=1, t 
IN= I*  5 

IF  (N  (IN)  . NL.  1)  GO  TO  15 
IS=I*  20 

S (I S)=X  (K)  +H  (3)  *T  (K) 

K = K*  1 

1 5 CONTINUE 
CrtLL  2HCK 
K 0 U N T=  K 0 U N T ♦ 1 
B (J)=S  (33) 

GO  TO  11 
C 

C NOW  CrlLCK  FOR  TERMINATION 

C 

14  IF  (KOUHT.GT.  N (14)  ) GO  TO  16 

C 

C CHECK  TO  SLt  IF  MINIMUM 

C HAS  OCCURSL  IN  PRESENT 

C DIRECTION 

C 

HL=H(2) 

I'U=ABS  (tl(J)-B  (2)  ) 

TL=  AHS  (B  (2)  -B  (1)  ) 

IF  (IL.GT.S(137)  ) GOTO  17 
i»M=  (11(2)  *H(1)  )/2.0 
K=  1 

DO  18  1=1,6 
1N=I*5 

IF  (N  (IN)  .NL.  1)  GO  TO  18 
IS=I*20 

S (IS)=X  (K)  *HM*T  (K) 
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K=K*  1 

15  CONTINUE 
CALL  ERCK 
KOUNT=KOUVT* 1 

E=  AbS  (5  (33) -ti  (2)  ) 

IF  (E.LE.S  (107)  ) GO  TO  20 
GO  TO  19 

17  IF  (TI).  GT.  S ( 1 37) ) GO  TO  19 
11M=  (ii  (3)  *H  (2)  )/2.  0 
K=1 

DO  21  1=1,  to 
IN= 1*5 

IF  (N  (IN)  . NE.  1)  GO  TO  21 
I S=  I ♦ 2 0 

S (IS)=X  (K)  *ri  M*T  (K) 

K = K ♦ 1 

21  CONTINUE 
CALL  ERCK 
KCUNT=KCUNT*  1 
E=  AbS  (S  ( 3 3)  - B (2)  ) 

If  (E. LI. S (107) ) GO  TO  20 

19  EELT=  DLLT/2. 

C 

C THE  MINIMUM  HAS  NOT  YET 

C BEEN  FOUND  IN  THIS  DIRECTION. 

C REDUCE  CLUSTER  SIZE  DY  HALF 

C S CONTINUE  SEARCH. 

C 

EX0  = H (2) 
b (1)  = EX0-DELT 
H (3)  =DX3*DbLT 
GO  TO  22 

20  CONTINUE 
C 

C THE  KEAN  SQUARED  ERROR  HAS  BEEN 

C SUFFICIENTLY  REDUCED  IN  THIS 

C DIRECTION.  IF  THIS  ERROR  IS 

C SMALL  ENOUGH  THE  SUBROUTINE 

C WILL  RETURN.  IP  NOT, 

C THE  SUBROUTINE  WILL  CONTINUE 

C SEARCH  IN  A NEW  DIRECTION. 

C 

IF  (S  (33)  .LE.  S-(93) ) RETURN 
K=  1 

DO  23  1=1,6 
IN=I*5 

IF  (N  (IN)  .NE.  1)  GO  TO  23 
X (K)=S  (1*20) 

K=K*  1 

23  CONTINUE 
GO  TO  24 
5 CONTINUE 
PRINT  25 

25  FORMAT  (111  , biiOUTR ANGE) 

STOP 

16  DO  26  1=1,6 
xS=I*20 
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3 . 4 SAMPLE  OUTPUT 

The  following  pages  provide  a sample  output  (as  run 
in  the  time  sharing  mode) . The  input  data  file  is  also 
provided,  although  the  data  file  structure  will  be  dis- 
cussed more  thoroughly  in  section  3.5.  Referring  to  the 
actual  output,  the  message  ERCK  is  printed  to  show  that 
the  program  is  actually  running  (as  opposed  to  being 
entangled  in  an  infinite  loop) . The  first  three  printed 
values  are  the  elevation  angles  for  minimum  group  delays 
at  the  first  three  sounding  frequencies.  Several  ERCK 
messages  are  then  printed  again  to  show  that  the  program 
is  still  running.  The  message  " *TIME  LIMIT*"  is  generated 
by  the  time  sharing  terminal  to  show  that  a specified  amount 
of  computer  time  has  been  used,  and  that  unless  additional 
time  is  allocated,  the  terminal  will  log  off.  The  command 
"T,20" , therefore,  is  entered  from  the  terminal  in  order 
to  allow  the  program  to  continue.  The  next  block  of  infor- 
mation is  a summary  of  the  results  computed  by  the  0th 
iteration.  The  message  " *C0NCAVE*"  indicates  that  certain 
necessary  conditions  for  the  concavity  of  the  MSE  surface 
have  been  satisfied.  The  present  ionostate  (which,  in  this 
case,  is  the  initial  ionostate)  is  then  printed,  followed 
by  a listing  of  the  components  of  the  unit  vector  along  the 
search  direction.  Since  there  are  three  variable  ionostate 
components  in  this  case  (GRB , GRM,  and  GFC) , there  are 
3 non-zero  components.  Finally,  the  value  of  the  mean 
squared  error  (MSE)  is  printed.  The  next  set  of  "ERCK" 
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messages  imply  that  the  program  is  looking  for  a minimum 
MSE  by  varying  the  variable  ionostate  components  in  a 
manner  consistent  with  the  previously  determined  search 
direction. 

Finally,  a local  minimum  is  found  and  the  elevation 
angles  are  for  the  minimum  group  paths  at  the  three 
sounding  frequencies  are  printed.  The  "*DS="  messages 
indicate  that  some  difficulty  was  experienced  in  computing 
the  derivatives  which  are  necessary  in  the  determination 
of  the  next  search  direction.  These  difficulties  arise 
when  the  finite  difference  method,  which  is  employed  to 
approximate  the  derivative,  finds  that  the  difference  in 
minimum  group  paths  is  about  the  same  size  as  the  uncer- 
tainty in  the  minimum  group  path  (which  we  believe  is 
on  the  order  of  a few  tenths  of  a kilometer  for  integration 
step  sizes  ranging  from  a few  hundred  meters  to  a few 
kilometers) . In  this  case  the  derivative  may  be  very 
unreliable  because  of  the  relatively  large  influence  of 
"noise".  On  the  other  hand,  a difference  in  minimum 
group  paths  which  is  many  orders  of  magnitude  larger  than 
this  uncertainty  is  objectionable  because  it  suggests  that 
the  finite  difference  approximation  may  be  too  crude. 
Therefore  occasional  adjustment  of  the  differentiation 
step  size,  as  seen  in  this  example,  is  performed  auto- 
matically. 

The  final  block  of  information  is  the  result  of  the 
first  complete  iteration.  It  can  be  seen  that  the  variable 
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ionostate  components  have  been  adjusted  and  a new  search 
direction  has  been  computed.  The  mean  squared  error  is  now 
much  smaller  and  the  final  message  "**CRASH005"  indicates 
that  the  program  was  terminated  because  the  number  of 
search  directions  was  specified  as  1 in  the  input  file. 

3.5  DATA  FILE  STRUCTURE 

Although  the  data  file  provided  in  section  3.4  suggests 
an  indication  of  some  typical  values  for  the  data  parameters , 
a more  detailed  discussion  may  be  beneficial  to  the  user 
of  this  program.  The  following  page,  therefore,  consists 
of  a computer  print  out  of  an  information  file  (’l)  which 
provides  an  overview  of  the  allocations  of  the  S and  N 
arrays.  These  two  arrays  serve  the  dual  purposes  of  pro- 
viding both  the  data  input  to  the  program  and  acting  as 
common  arrays  to  be  passed  between  the  various  subroutines. 
Accordingly,  the  array  elements  which  are  used  as  data 
input  are  denoted  by  the  usage  symbol  E (external)  and  the 
elements  used  for  transferring  information  between  sub- 
routines are  denoted  by  the  symbol  I (internal) . 

In  the  following  discussion,  however,  we  will  be  con- 
cerned only  with  the  external  elements  of  the  S and  N arrays. 
In  particular,  it  is  our  intention  to  mention  some  of  the 
specific  considerations  which  must  be  made  in  the  assignment 
of  values  to  these  elements.  Elements  S(2)-S(7),  for 
example,  are  self-explanatory  although  the  user  should  be 
aware  that  when  less  than  the  maximum  number  of  sounding 
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frequencies  (6)  are  used,  the  values  should  be  entered  in 
consecutive  elements,  starting  from  S(2).  The  next  group 
of  elements,  S(8)-S(13),  are  minimum  group  paths,  and  while 
these  elements  are  entered  in  the  same  order  as  S(2)-S(7), 
they  are  neither  required  nor  used  in  the  synthesis  mode  of 
the  program.  S(14)-S(19)  are  then  the  corresponding  eleva- 
tion angles  which  produce  the  minimum  group  paths  in  S(8)-S(13). 
In  general  these  angles  are  unknown,  but  are  used  only  as 
starting  values  for  subroutine  STAR.  It  has  been  our  practice 
to  use  the  synthesis  mode  to  obtain  this  information,  but 


any  reasonable  "guesses"  should  be  sufficient.  The  Runge- 
Kutta  step  size,  S(20),  however,  is  rather  critical  since 
an  excessively  large  value  will  cause  significant  errors 
in  the  ray  tracing  subroutine,  RAYT,  while  an  excessively 
small  value  will  cause  the  expenditure  of  substantial 
amounts  of  computer  time.  Typical  values  for  this  parameter 
would  be  in  the  approximate  range  of  . 1 to  5 km.  The  iono- 
state  is  entered  in  the  order:  RB,  RM,  FC,  GRB , GRM,  GFC. 

In  the  inversion  mode,  this  ionostate  corresponds  to  the 
initial  guess  of  the  correct  ionostate,  while  in  the  synthesis 
mode  these  parameters  describe  the  ionosphere  in  which  the 
backscatter  leading  edge  is  simulated.  The  elements  S(83) 
and  S(84)  are  used  in  subroutine  STAR.  S(83)  determines 
when  the  minimum  group  path  has  been  found  (see  discussion 
of  subroutine  STAR).  For  obvious  reasons,  S(83),  must  not 
be  too  much  smaller  than  the  uncertainty  in  ray  tracing. 

The  initial  star  size,  S(84),  was  also  discussed  in  connection 
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with  subroutine  STAR.  S(87)-S(92)  are  differentiation 
step  sizes  which  are  used  in  approximating  the  gradient 
of  the  MSE  by  the  finite  difference  method.  Only  those 
elements  corresponding  to  variable  ionostate  components 
need  to  be  specified.  It  should  also  be  mentioned  that 
these  elements  may  be  adjusted  by  the  program  if,  based 
upon  certain  criteria,  it  would  appear  that  the  derivatives 
may  be  unreliable,  as  discussed  previously  in  connection 
with  the  sample  output.  S(93)  provides  one  of  the  termina- 
tion criteria  for  subroutine  STEEP  and  its  value  is  somewhat 
subjective.  Our  typical  value  of  1 km£  is  sufficient  be- 
cause we  have  used  simulated  data  in  an  LQP  model  ionosphere 
and  thus  have  known  that  some  ionostate,  at  least  in  theory, 
would  result  in  an  MSE  of  0 km2.  In  a more  realistic  case, 
a final  MSE  of  1 km2  or  less  may  not  be  attainable  and 
thus  a larger  value  for  S(93)  may  be  used.  However,  it 
should  be  emphasized  that  even  if  S(93)  is  too  small, 
other  parameters  will  eventually  cause  program  termination 
and  therefore  there  is  no  real  danger  of  performing  an 
infinite  search  for  an  unrealizable  MSE.  S(95)-S(106) 
define  a finite  domain  in  ionostate  space  over  which  the 
search  for  the  optimum  ionostate  is  to  be  performed. 

Again,  only  the  bounds  on  variable  ionostate  components 
need  to  be  specified.  It  is,  of  course,  necessary  that 
the  initial  ionostate  be  contained  in  this  domain.  In 
the  event  that  the  program  begins  to  use  ionostates  outside 
of  this  domain  a termination  will  occur  which  usually  implies 
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3.6  CRASH  CODES 

The  following  codes  are  used  to  indicate  the  reason 
for  program  termination 

TABLE  2:  CRASH  CODE  SUMMARY 


CODE 


IMPLICATION 


NSTP>N (1) 

Too  many  iterations  in  RAYT. 

Try  larger  value  of  N(l) 

CRASH001 

Too  many  calls  to  RAYT. 

Try  large  value  for  N(5) 

CRASH002 

No  minimum  group  path  could  be 
found  for  the  given  sounding 
frequency.  Check  data  input 
file  for  mistakes.  If  file  is 
correct,  different  sounding 
frequencies  may  b':  necessary. 

CRASH004 

No  differentiation  step  size 
specified.  Check  data  input 
file  to  make  sure  that  for 
every  variable  ionostate  com- 
ponent there  is  a differentia- 
tion step  size  in  S(87)-S(92) 

CRASH005 

Program  terminated  upon  reaching 
specified  number  of  search 
directions  without  finding 
specified  minimum  mean  squared 
error. 

OUTRANGE 

I 

Ionostate  is  outside  of  speci- 
fied domain.  This  domain  may 
have  to  be  made  larger,  and/or 
a different  starting  ionostate 
may  be  required. 

ITERATION 

Too  many  iterations  in  subrou- 
tine STEEP. 

Try  increasing  N(14). 
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4.  CONCLUSION 

As  a result  of  our  emphasis  on  applications,  examples, 
and  documentation  in  the  previous  sections,  there  are  a 
few  considerations  of  a general  nature  which  still  remain 
to  be  discussed.  In  addition,  we  would  like  to  make  some 
specific  recommendations  for  possible  improvements  in  the 
program. 

The  most  persistent  problem  which  we  have  considered 
is  the  problem  of  uniqueness  in  computing  an  optimum  iono- 
state.  The  question  of  whether  or  not  two  different 
ionostates  could  result  in  identical  backscatter  leading 
edges  is  still  unresolved.  On  a qualitative  basis,  we  have 
discussed  (see  section  2.4)  the  possibility  of  various 
ionostate  components  having  a mutually  offsetting  effect 
and  it  would,  at  least,  therefore  seem  plausible  that  there 
is  not  a one-to-one  mapping  between  the  backscatter  leading 
edge  and  a corresponding  ionostate.  However,  even  if  there 
were  such  a one-to-one  mapping  it  could  still  be  possible 
that  the  mean  squared  error  hypersurface  would  have  more 
than  one  local  minimum  as  a function  of  ionostate.  In 
this  case,  the  values  at  each  of  the  local  minima  would 
be  different  and  the  objective,  of  course,  would  be  to 
find  the  local  minimum  corresponding  to  the  least  mean 
squared  error.  Using  the  present  form  of  this  program  in 
order  to  find  the  smallest  local  minimum  would  then  require 
running  the  inversion  mode  for  many  different  starting 


97 


values  for  the  initial  ionostate.  This  problem  (i.e. 
sensitivity  to  starting  values)  is  one  that  is  common 
to  all  of  the  standard  minimization  procedures  that  we 
are  presently  acquainted  with.  We  believe,  however,  that 
the  choice  of  a starting  value  would  be  facilitated  by  the 
generation  of  a library  of  the  backscatter  leading  edges 
produced  by  ray  tracing  in  LQP  model  ionospheres  based 
upon  various  ionostates.  The  selection  of  a particular 
starting  value  for  the  ionostate  may  also  be  facilitated 
through  the  use  of  the  N(21)  flag  (see  section  3.5).  We 
have  used  this  procedure  to  our  advantage  in  obtaining  a 
number  of  the  results  presented  in  section  2.  In  particu- 
lar, setting  N(21)=l  and  using  various  ionostates  in  S(21)- 
S(26)  can  provide  a rough  picture  of  the  MSE  hypersurface. 

Another  problem,  which  was  previously  discussed,  is 
the  selection  of  the  Runge-Kutta  step  size.  The  present 
ray  tracing  subroutine  (RAYT)  uses  a fixed  value,  S(20), 
for  the  step  size.  A more  efficient  procedure  would  be 
to  approximate  both  the  truncation  and  round-off  errors 
in  each  integration  step,  and  adjust  the  step  size  accord- 
ingly. There  is,  however,  some  utilization  of  this  idea 
in  the  present  program  through  the  incorporation  of  separate 
subroutines  for  free  space  and  ionospheric  ray  tracing. 

In  spite  of  these  problems,  a number  of  successful 


computer  runs  have  confirmed  the  feasibility  of  this  pro- 
cedure, and  resulted  in  ionostates  which  were  correct  to 
within  the  tolerances  specified  in  the  data  input  file. 
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